function IAR=IntActReg(Y,X,M,AT,F,NL,varargin);%


K=size(M,2)/size(X,2);%Number of interaction terms
[T,N]=size(X);


%for k=1:K
 %XX(:,(k-1)*N+1:k*N)= M(:,(k-1)*N+1:k*N).*X;
%end;
%size(XX)
XX=kron(ones(1,K),X); 
size(XX);


if AT>0
    XX=[M.*XX,M(:,N*(K-AT)+1:size(M,2))];
    %size(XX)
    %XX=[M.*XX,NL];
    XX=reshape(XX,T*N,K+AT);
else
    XX=M.*XX;
    XX=reshape(XX,T*N,K);
    size(XX);
end

if F==1
FE=kron(eye(N),ones(T,1));FE=FE(:,1:N-1);
XX=[ones(size(XX,1),1),XX,FE];
%Reg=NewOLS(Y,[XX,ones(size(XX,1),1),FE]);

elseif F==2
    FE=kron(ones(N,1),eye(T));FE=FE(:,1:T-1);
    XX=NewOLS(Y,[ones(size(XX,1),1),XX,FE]);
    %Reg=NewOLS(Y,[XX,ones(size(XX,1),1),FE]);
elseif F==3
FE1=kron(eye(N),ones(T,1));FE1=FE1(:,1:N-1);
FE2=kron(ones(N,1),eye(T));FE2=FE2(:,1:T-1);
XX=[ones(size(XX,1),1),XX,FE1,FE2];
%Reg=NewOLS(Y,[ones(size(XX,1),1),XX,FE1,FE2]);
else
XX=[ones(size(XX,1),1),XX];    
%Reg=NewOLS(Y,[ones(size(XX,1),1),XX]);

end;

Reg=NewOLS(Y,[XX]);
IAR=Reg;

if nargin==7;
        optionalarg=varargin{1};
  switch optionalarg
    case 'cluster'
     clusters=kron((1:N)',ones(T,1));    
     for n=1:size(Y,2);
     ClusterReg{n}=clusterreg(Y(:,n),[XX],clusters);
     end;
     IAR.ClusterReg=ClusterReg;
      case 'twocluster'
          cluster_regions=kron((1:N)',ones(T,1));
          cluster_time=kron(ones(N,1),(1:T)');
          for n=1:size(Y,2);
          ClusterReg{n}=clusterreg(Y(:,n),[XX],cluster_regions,cluster_time);
          end;
          IAR.ClusterReg=ClusterReg;    
          
    otherwise
     disp('invalid option')
  end;

end;



Fitted_b=zeros(size(Y));
for k=1:K
    for l=1:size(Y,2);
Fitted_b(:,l)=Fitted_b(:,l)+IAR.b(k,l)*vec2mat(M(:,(k-1)*N+1:k*N),1);
end;%l
end;%k

IAR.Fitted_b=Fitted_b;






